Integrate Virus scRNA-seq data

library(Seurat)
library(ggplot2)
library(dplyr)
library(harmony)
library(clusterProfiler)
library(enrichplot)
library(msigdbr)
library(stringr)
library(tidyr)
library(patchwork)
library(Azimuth)
set.seed(2025)

load data

out_path <- "~/projects/2024/RA/Flu_integration/"
# load E_MTAB_12392
in_path <- "~/projects/2024/RA/E_MTAB_12392/out/"
E_MTAB_12392_scdata <- readRDS(paste0(in_path,"E_MTAB_12392_scdata.rds"))
E_MTAB_12392_metadata <- readRDS(paste0(in_path,"E_MTAB_12392_metadata.rds"))
length(E_MTAB_12392_scdata)

[1] 18

# load GSE164381
in_path <- "~/projects/2024/RA/GSE164381_RAW/out/"
GSE164381_scdata <- readRDS(paste0(in_path,"GSE164381_scdata.rds"))
GSE164381_metadata <- readRDS(paste0(in_path,"GSE164381_metadata.rds"))
length(GSE164381_scdata)

[1] 11

select passed cutoff cells in each cohort

# sample name
d_name_1 <- unique(E_MTAB_12392_metadata$orig.ident)
d_name_2 <- unique(GSE164381_metadata$orig.ident)
d_name <- c(d_name_1,d_name_2)
# cell list
cell_list_1 <- rownames(E_MTAB_12392_metadata)
cell_list_2 <- rownames(GSE164381_metadata)
head(cell_list_1)

[1] "S1_AAACCTGAGATAGTCA-1" "S1_AAACCTGAGGAACTGC-1" "S1_AAACCTGGTGTCGCTG-1"
[4] "S1_AAACGGGTCACATGCA-1" "S1_AAACGGGTCCAACCAA-1" "S1_AAACGGGTCCAATGGT-1"

head(cell_list_2)

[1] "GSM5008746_AAACCTGAGTGACATA-1" "GSM5008746_AAACCTGCATTAGCCA-1"
[3] "GSM5008746_AAACCTGGTCTCGTTC-1" "GSM5008746_AAACCTGGTGCCTGGT-1"
[5] "GSM5008746_AAACCTGTCACTTACT-1" "GSM5008746_AAACGGGAGACAGACC-1"

cell_list <- c(cell_list_1, cell_list_2)
# merge
scdata <- c(E_MTAB_12392_scdata, GSE164381_scdata)
Sobj <- merge(scdata[[1]], y = scdata[-1], add.cell.ids = d_name)
Sobj[["percent.mt"]] <- PercentageFeatureSet(Sobj, pattern = "^MT-")
# select passed cells 
Sobj <- subset(Sobj, cells = cell_list)
length(cell_list)

[1] 39055

dim(Sobj)

[1] 36601 39055

normalize, scale, PCA

Sobj <- NormalizeData(Sobj, normalization.method = "LogNormalize", scale.factor = 10000)

Normalizing layer: counts.S1
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S3
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S11
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S13
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S14
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S16
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S17
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S18
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S19
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S20
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S21
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S22
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S23
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S24
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S25
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S26
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S27
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S28
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008746
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008747
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008748
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008749
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008750
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008751
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008752
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008753
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008754
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008755
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008756
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Sobj <- FindVariableFeatures(Sobj, selection.method = "vst", nfatures = 2000)

Finding variable features for layer counts.S1
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S3
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S11
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S13
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  at  -0.31854
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  radius  0.00030656
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  all data on boundary of neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -0.31854
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.017509
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  zero-width neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.090619
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S14
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S16
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S17
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S18
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S19
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S20
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -1.9345
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.30103
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  6.7402e-15
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S21
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -1.752
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.30103
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1.0124e-14
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S22
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  at  -0.91849
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  radius  0.00023715
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  all data on boundary of neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -0.91849
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.0154
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  zero-width neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.031008
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S23
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S24
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -2.1206
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.30103
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  8.0208e-15
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S25
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -1.2225
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.31945
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  7.6188e-28
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.031008
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S26
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  at  -0.62058
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  radius  0.00034289
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  all data on boundary of neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -0.62058
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.018517
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  zero-width neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.031008
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S27
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  at  -1.0144
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  radius  0.0002087
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  all data on boundary of neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -1.0144
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.014447
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  zero-width neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.031008
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S28
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -2.2407
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.49766
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1.8412e-14
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.090619
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008746
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008747
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008748
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008749
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008750
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008751
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008752
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008753
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008754
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008755
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008756
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Sobj <- ScaleData(Sobj)

Centering and scaling data matrix

  |                                                                                             
  |                                                                                       |   0%
  |                                                                                             
  |============================================                                           |  50%
  |                                                                                             
  |=======================================================================================| 100%

Sobj <- RunPCA(Sobj)

PC_ 1 
Positive:  MT-CO2, MT-CO3, MT-CYB, CD83, MALAT1, MS4A1, MT-CO1, CD74, MT-ATP6, SNX9 
       BACH2, HLA-DQB1, ZNF331, SLC2A3, CD79A, KLF6, HLA-DRA, FCER2, MT-ND5, CD69 
       PHACTR1, HLA-DRB1, LY9, IGHM, NEIL1, IL4R, SERPINB9, CCR7, HLA-DQA1, IGHD 
Negative:  GAPDH, ACTG1, S100A4, PFN1, SRGN, ANXA1, ACTB, S100A6, AIF1, CLIC1 
       HCST, TPI1, CTSW, SH3BGRL3, TXN, TALDO1, TYROBP, CFL1, GSTP1, IFITM2 
       ENO1, IGFBP7, FCER1G, VIM, S100A10, MYL12A, PRSS57, CD99, PPIB, ARPC1B 
PC_ 2 
Positive:  HLA-DRA, CD74, HLA-DPA1, HLA-DPB1, RPS8, RPLP0, HLA-DRB1, RPLP1, EEF1B2, RACK1 
       RPS18, RPS4X, RPL5, RPS12, RPS6, RPS2, EEF1A1, RPL7A, RPL10, RPS3A 
       CD79A, RPS14, RPS24, FTL, RPS3, HLA-DQB1, SNHG29, RPL7, HLA-DRB5, FTH1 
Negative:  GNLY, NKG7, CST7, GZMB, CCL5, GZMA, CTSW, CD247, CD7, GZMM 
       FGFBP2, PRF1, KLRD1, GZMH, SPON2, KLRB1, TYROBP, CLIC3, CCL4, HCST 
       ZAP70, FCGR3A, IL2RB, IFITM1, IL32, S1PR5, KLRF1, TRBC1, PRKCH, MATK 
PC_ 3 
Positive:  PRSS57, H1F0, CYTL1, SMIM24, EGFL7, GATA2, AL157895.1, BEX3, GIHCG, NFE2 
       SLC40A1, CPA3, SOX4, MT-CO1, MT-CO3, CD34, EREG, EMID1, MT-CYB, FCER1A 
       MPP1, CRHBP, LAPTM4B, CDK6, IL18, ERG, DEPTOR, LINC02573, NPR3, MALAT1 
Negative:  TMSB4X, HLA-DPA1, TMSB10, B2M, CD79A, ACTB, RPS18, LTB, FTL, RPL10 
       HLA-DPB1, CYBA, EEF1A1, RPLP1, RPS3, CD74, RPS12, HLA-DRA, LSP1, CORO1A 
       RPL7A, HLA-DRB1, CFL1, JUNB, RPS4X, PFN1, RPS14, RPS6, EEF1B2, DUSP1 
PC_ 4 
Positive:  RPS3, EEF1A1, RPS12, RPL10, RPL7A, RPS3A, RPS4X, RPS14, RPL5, RPS8 
       RACK1, RPS6, RPS18, EEF1B2, RPS2, RPLP0, TPT1, RPLP1, RPS24, CTSW 
       NKG7, RPL7, CST7, NPM1, SNHG29, GAS5, CD7, GNLY, HNRNPA1, GZMA 
Negative:  LYZ, CST3, IFI30, S100A9, CLEC7A, FCN1, CPVL, SERPINA1, S100A8, CSTA 
       PLAUR, LGALS3, CFP, LGALS2, CFD, S100A11, TNFAIP2, CD68, G0S2, MAFB 
       AIF1, S100A12, TGFBI, MS4A6A, VCAN, TYMP, SHTN1, AC020656.1, HMOX1, PLBD1 
PC_ 5 
Positive:  TXNDC5, JCHAIN, TNFRSF17, MZB1, IGHA1, CD27, AQP3, PDIA4, DERL3, IGHG1 
       IGHA2, ITM2C, SEC11C, MTRNR2L12, FKBP11, IGKC, IGHG2, PRDM1, SRM, CHPF 
       CLPTM1L, MAN1A1, SDF2L1, CD38, MANF, LMAN1, PRDX4, STT3A, SSR3, GPRC5D 
Negative:  FTH1, LYZ, CST3, RPS14, RPS3A, RPL10, S100A9, FCN1, RPS12, RPS8 
       RPS2, EEF1A1, TYROBP, AIF1, TPT1, RPS24, FCER1G, SERPINA1, S100A8, CSTA 
       CLEC7A, FTL, RPLP1, RPL7A, RPS4X, RPS3, RPL5, PLAUR, LST1, CPVL 

select cutoff of dimension.

ElbowPlot(Sobj)

PCA space

DimPlot(Sobj, group.by = "orig.ident")

UMAP

Sobj <- RunUMAP(Sobj, dims = 1:15, reduction = "pca", reduction.name = "umap.unintegrated")

20:06:19 UMAP embedding parameters a = 0.9922 b = 1.112
20:06:19 Read 39055 rows and found 15 numeric columns
20:06:19 Using Annoy for neighbor search, n_neighbors = 30
20:06:19 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:06:24 Writing NN index file to temp file /tmp/RtmpBxlMZd/file5c7036bb70824
20:06:24 Searching Annoy index using 1 thread, search_k = 3000
20:06:38 Annoy recall = 100%
20:06:39 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:06:43 Initializing from normalized Laplacian + noise (using RSpectra)
20:06:46 Commencing optimization for 200 epochs, with 1667658 positive edges
20:06:46 Using rng type: pcg
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:07:08 Optimization finished

p0 <- DimPlot(Sobj, reduction = "umap.unintegrated", group.by = c("orig.ident", "Batch"))

Warning: The following requested variables were not found: Batch

p0

Harmony batch correction

Sobj <- RunHarmony(Sobj,"orig.ident", plot_convergence = T, theta = 3)

Transposing data matrix
Initializing state using k-means centroids initialization
Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 5 iterations

harmony_embeddings <- Embeddings(Sobj, "harmony")
harmony_embeddings[1:5, 1:5]

                       harmony_1 harmony_2  harmony_3   harmony_4   harmony_5
S1_AAACCTGAGATAGTCA-1   4.701568 -1.708613  1.3231752  -2.8416889  -0.1320850
S1_AAACCTGAGGAACTGC-1   2.451071  1.390329 -1.8605446   1.6223532   0.2989081
S1_AAACCTGGTGTCGCTG-1   4.284953  2.023095 -0.8945279   1.5819618  -1.1558085
S1_AAACGGGTCACATGCA-1   3.936914  1.105024 -1.2764343   0.8342876  -1.0429716
S1_AAACGGGTCCAACCAA-1 -23.461463  1.001002 -3.2267307 -25.5267511 -14.3578007

DimPlot(Sobj, reduction = "harmony", group.by = "orig.ident")

UMAP with harmony

Sobj <- RunUMAP(Sobj, dims = 1:15, reduction = "harmony")

20:07:59 UMAP embedding parameters a = 0.9922 b = 1.112
20:07:59 Read 39055 rows and found 15 numeric columns
20:07:59 Using Annoy for neighbor search, n_neighbors = 30
20:07:59 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:08:04 Writing NN index file to temp file /tmp/RtmpBxlMZd/file5c7034015e3f0
20:08:04 Searching Annoy index using 1 thread, search_k = 3000
20:08:19 Annoy recall = 100%
20:08:20 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:08:23 Initializing from normalized Laplacian + noise (using RSpectra)
20:08:26 Commencing optimization for 200 epochs, with 1710788 positive edges
20:08:26 Using rng type: pcg
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:08:48 Optimization finished

DimPlot(Sobj, reduction = "umap", group.by = "orig.ident")

Check whether the harmony batch correction work well or not.

Idents(Sobj) <- Sobj$orig.ident
rename_map <- c(
    setNames(rep("E_MTAB_12392", length(d_name_1)), d_name_1),
    setNames(rep("GSE164381", length(d_name_2)), d_name_2)
)
Sobj <- RenameIdents(Sobj, rename_map)
unique(Idents(Sobj))

[1] E_MTAB_12392 GSE164381   
Levels: E_MTAB_12392 GSE164381

Sobj$batch <- Idents(Sobj)

DimPlot(Sobj, group.by = "batch")

DimPlot(Sobj, reduction = "umap", group.by = "batch")

Sobj

An object of class Seurat 
36607 features across 39055 samples within 2 assays 
Active assay: RNA (36601 features, 2000 variable features)
 59 layers present: counts.S1, counts.S3, counts.S11, counts.S13, counts.S14, counts.S16, counts.S17, counts.S18, counts.S19, counts.S20, counts.S21, counts.S22, counts.S23, counts.S24, counts.S25, counts.S26, counts.S27, counts.S28, counts.GSM5008746, counts.GSM5008747, counts.GSM5008748, counts.GSM5008749, counts.GSM5008750, counts.GSM5008751, counts.GSM5008752, counts.GSM5008753, counts.GSM5008754, counts.GSM5008755, counts.GSM5008756, data.S1, data.S3, data.S11, data.S13, data.S14, data.S16, data.S17, data.S18, data.S19, data.S20, data.S21, data.S22, data.S23, data.S24, data.S25, data.S26, data.S27, data.S28, data.GSM5008746, data.GSM5008747, data.GSM5008748, data.GSM5008749, data.GSM5008750, data.GSM5008751, data.GSM5008752, data.GSM5008753, data.GSM5008754, data.GSM5008755, data.GSM5008756, scale.data
 1 other assay present: HTO
 4 dimensional reductions calculated: pca, umap.unintegrated, harmony, umap

Sobj <- FindNeighbors(Sobj, dims = 1:15, reduction = "harmony")

Computing nearest neighbor graph
Computing SNN

Sobj <- FindClusters(Sobj, reduction = "umap", resolution = 1)

Warning: The following arguments are not used: reduction
Warning: The following arguments are not used: reduction

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 39055
Number of edges: 1156538

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8409
Number of communities: 18
Elapsed time: 10 seconds

Sobj@reductions

$pca
A dimensional reduction object with key PC_ 
 Number of dimensions: 50 
 Number of cells: 39055 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

$umap.unintegrated
A dimensional reduction object with key umapunintegrated_ 
 Number of dimensions: 2 
 Number of cells: 39055 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

$harmony
A dimensional reduction object with key harmony_ 
 Number of dimensions: 50 
 Number of cells: 39055 
 Projected dimensional reduction calculated:  TRUE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

$umap
A dimensional reduction object with key umap_ 
 Number of dimensions: 2 
 Number of cells: 39055 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

DimPlot(Sobj, reduction = "umap", group.by = "RNA_snn_res.1")

Sobj <- JoinLayers(Sobj)

DimPlot(Sobj, reduction = "umap", label = T)

load reference data

SLE_Sobj <- readRDS("~/projects/2024/RA/Integration_SLE_Sobj.rds")
SLE_Sobj <- RunUMAP(SLE_Sobj, reduction = "harmony",dims = 1:15, return.model = TRUE)

UMAP will return its model
20:10:08 UMAP embedding parameters a = 0.9922 b = 1.112
20:10:08 Read 133389 rows and found 15 numeric columns
20:10:08 Using Annoy for neighbor search, n_neighbors = 30
20:10:08 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:10:24 Writing NN index file to temp file /tmp/RtmpBxlMZd/file5c703562f00f7
20:10:24 Searching Annoy index using 1 thread, search_k = 3000
20:11:27 Annoy recall = 100%
20:11:29 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:11:34 Initializing from normalized Laplacian + noise (using RSpectra)
20:11:47 Commencing optimization for 200 epochs, with 5832214 positive edges
20:11:47 Using rng type: pcg
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:13:09 Optimization finished

DimPlot(SLE_Sobj, reduction = "umap", label = T, repel = T)

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

DimPlot(SLE_Sobj, reduction = "umap.unintegrated")

Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

# if (!file.exists(paste0(out_path,"Integration_Virus.RData"))) {
#   save.image(file = paste0(out_path,"Integration_Virus.RData"))
# }

calculate transfer anchors

anchors <- FindTransferAnchors(reference = SLE_Sobj,
                               query = Sobj,
                               normalization.method = "LogNormalize",
                               reference.reduction = "harmony",
                               dims = 1:15)

Projecting cell embeddings
Finding neighborhoods
Finding anchors
    Found 6697 anchors

Mapping

Sobj <- MapQuery(anchorset = anchors,
                 query = Sobj,
                 reference = SLE_Sobj,
                 refdata = list(celltype = "anno_text"),
                 reference.reduction = "harmony",
                 reduction.model = "umap")

Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Predicting cell labels
Warning: Layer counts isn't present in the assay object; returning NULL
Warning: Layer counts isn't present in the assay object; returning NULL
Warning: Layer counts isn't present in the assay object; returning NULL

  |                                                  | 0 % ~calculating  

Integrating dataset 2 with reference dataset
Finding integration vectors
Integrating data

  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  

Computing nearest neighbors
Running UMAP projection
20:16:20 Read 39055 rows
20:16:20 Processing block 1 of 1
20:16:20 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:16:21 Initializing by weighted average of neighbor coordinates using 1 thread
20:16:21 Commencing optimization for 67 epochs, with 1171650 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:16:30 Finished

Idents(Sobj) <- Sobj$predicted.celltype
Idents(Sobj) <- factor(x = Idents(Sobj), levels = c("Naive B", "Activated B", "Unswitched memory B", "CD27+ memory B", "Plasmablast", "Plasma", "T cell like B", "Myeloid like B", "Megakaryocyte like B"))
levels(Idents(Sobj))

[1] "Naive B"              "Activated B"          "Unswitched memory B"  "CD27+ memory B"      
[5] "Plasmablast"          "Plasma"               "T cell like B"        "Myeloid like B"      
[9] "Megakaryocyte like B"

levels(SLE_Sobj$anno_text)

 [1] "Naive B"              "Activated B"          "Unswitched memory B"  "CD27+ memory B"      
 [5] "Atypical memory B"    "INF-induced memory B" "Plasmablast"          "Plasma"              
 [9] "T cell like B"        "Myeloid like B"       "Megakaryocyte like B"

Sobj$predicted.celltype <- Idents(Sobj)

table(Sobj$predicted.celltype)

             Naive B          Activated B  Unswitched memory B       CD27+ memory B 
               29127                 2631                   53                 1264 
         Plasmablast               Plasma        T cell like B       Myeloid like B 
                  21                  387                 3885                 1684 
Megakaryocyte like B 
                   3 

Idents(Sobj) <- Sobj$predicted.celltype
Sobj <- subset(Sobj, idents = setdiff(levels(Sobj), "Megakaryocyte like B"))
table(Sobj$predicted.celltype)

            Naive B         Activated B Unswitched memory B      CD27+ memory B 
              29127                2631                  53                1264 
        Plasmablast              Plasma       T cell like B      Myeloid like B 
                 21                 387                3885                1684 

Idents(Sobj) <- Sobj$RNA_snn_res.1
# p1 <- DimPlot(Sobj, reduction = "umap", group.by = "RNA_snn_res.1")
p2 <- DimPlot(Sobj, reduction = "umap", label = T)
Idents(Sobj) <- Sobj$predicted.celltype
p3 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "ref.umap") + NoLegend()
p4 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "ref.umap")
p5 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "umap") + NoLegend()
p6 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "umap")

E_MTAB_12392_metadata$anno_ABC_LRB <- NA

E_MTAB_12392_metadata$anno_ABC_LRB[E_MTAB_12392_metadata$LSP1_score_pos==T & E_MTAB_12392_metadata$ABC_score_pos==F] <- "LRB"
E_MTAB_12392_metadata$anno_ABC_LRB[E_MTAB_12392_metadata$LSP1_score_pos==F & E_MTAB_12392_metadata$ABC_score_pos==T] <- "ABC"
E_MTAB_12392_metadata$anno_ABC_LRB[E_MTAB_12392_metadata$LSP1_score_pos==T & E_MTAB_12392_metadata$ABC_score_pos==T] <- "LRB & ABC"
E_MTAB_12392_metadata$anno_ABC_LRB[E_MTAB_12392_metadata$LSP1_score_pos==F & E_MTAB_12392_metadata$ABC_score_pos==F] <- "Classical B"

unique(E_MTAB_12392_metadata$anno_ABC_LRB)

[1] "Classical B" "ABC"         "LRB & ABC"   "LRB"        

unique(GSE164381_metadata$anno_ABC_LRB)

[1] ABC         Classical B LRB         LRB & ABC  
Levels: LRB ABC LRB & ABC Classical B

anno_vars <- c("anno_ABC_LRB", "LSP1_score", "ABC_score", "LSP1_score_pos", "ABC_score_pos")

for (var in anno_vars) {
  l1 <- E_MTAB_12392_metadata[[var]]
  l2 <- GSE164381_metadata[[var]]
  
  if (is.factor(l1)) l1 <- as.character(l1)
  if (is.factor(l2)) l2 <- as.character(l2)
  
  names(l1) <- rownames(E_MTAB_12392_metadata)
  names(l2) <- rownames(GSE164381_metadata)
  
  combined_anno <- c(l1, l2)
  
  Sobj[[var]] <- combined_anno[colnames(Sobj)]
}
rm(l1, l2, combined_anno)

my_colors <- c("LRB" = "#ff0000", "ABC" = "#ffff00", "LRB & ABC" = "#00b0ff", "Classical B" = "#BFBFBF")
Idents(Sobj) <- Sobj$anno_ABC_LRB
p7 <- DimPlot(Sobj, reduction = "umap", cols = my_colors) + NoLegend()
DimPlot(Sobj, reduction = "umap", cols = my_colors)

celltype.prop <- as.data.frame(table(Sobj$predicted.celltype, Sobj$anno_ABC_LRB))
# head(celltype.prop)
celltype.prop <- celltype.prop %>%
    rename(celltype = Var1, `ABC LRB` = Var2) %>%
    group_by(celltype) %>%
    mutate(Percent = Freq / sum(Freq)*100)
ggplot(celltype.prop, aes(x = celltype, y = Percent, fill = `ABC LRB`))+
    geom_bar(stat = "identity", width = 0.9) + 
    scale_fill_manual(values = my_colors) +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1), axis.title.x = element_blank(), plot.margin = margin(10, 10, 20, 10))

celltype_percent_wide <- celltype.prop %>%
  select(celltype, `ABC LRB`, Percent) %>%
  pivot_wider(names_from = `ABC LRB`, values_from = Percent)
head(celltype_percent_wide)

celltype_freq_wide <- celltype.prop %>%
  select(celltype, `ABC LRB`, Freq) %>%
  pivot_wider(names_from = `ABC LRB`, values_from = Freq)
head(celltype_freq_wide)

write.csv(celltype_percent_wide, paste0(out_path,"out/virus_barplot_percent.csv"))
write.csv(celltype_freq_wide, paste0(out_path,"out/virus_barplot_freq.csv"))

p2

p3

p4

p5

p6

p7
pdf(paste0(out_path, "Azimuth_predicted_umap_063025.pdf"), width = 4, height = 4)

p2
p3

p4
p5

p6
p7

dev.off()

png 
  2 

Idents(Sobj) <- Sobj$predicted.celltype
p6 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "umap")
pdf(paste0(out_path, "new_virus_umap_080425.pdf"), width = 4, height = 4)
p6
p7

dev.off()

null device 
          1 

p6_1 <- DimPlot(Sobj, label = FALSE, repel = TRUE, reduction = "umap") +
  labs(x = "", y = "") +
  theme(
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(),
    axis.ticks.y = element_line()
    )
tiff(filename = paste0(out_path,"umap_virus.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(p6_1)

dev.off()

null device 
          1 

p7_1 <- p7 +
  labs(x = "", y = "") +
  theme(
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(),
    axis.ticks.y = element_line()
    )
tiff(filename = paste0(out_path,"umap_LRB_ABC.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(p7_1)

dev.off()

null device 
          1 

table(Sobj$anno_ABC_LRB)

        ABC Classical B         LRB   LRB & ABC 
       8486       18588        6126        5852 

library(ggpointdensity)
library(scales)

df <- FetchData(Sobj, vars = c("umap_1", "umap_2", "LSP1_score_pos", "ABC_score_pos", "anno_ABC_LRB"))
df_LSP1 <- df %>% filter(LSP1_score_pos)
df_ABC <- df %>% filter(ABC_score_pos)

t1 <- ggplot(df_LSP1, aes(x = umap_1, y = umap_2)) +
  geom_pointdensity(method = "neighbors", size = 0.1) +
  scale_color_gradientn(
    colors = c("grey", "grey", "red"),
    values = rescale(c(0, 350, 1200)),
    limits = c(0, 1000)
  ) +
    theme_classic()

t2 <- ggplot(df_ABC, aes(x = umap_1, y = umap_2)) +
  geom_pointdensity(method = "neighbors", size = 0.1) +
  scale_color_gradientn(
    colors = c("grey", "grey", "red"),
    values = rescale(c(0, 350, 1200)),
    limits = c(0, 1000)
  ) +
    theme_classic()
t1

t2

pdf(paste0(out_path,"LRB_ABC_umap_density_080425.pdf"), width = 6, height = 5)
t1
t2

dev.off()

null device 
          1 

t1_1 <- t1 +
labs(x = "", y = "") +
theme(
  legend.position = "none",
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks.x = element_line(),
  axis.ticks.y = element_line()
  )
t2_1 <- t2 +
labs(x = "", y = "") +
theme(
  legend.position = "none",
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks.x = element_line(),
  axis.ticks.y = element_line()
  )
t1_1

t2_1

tiff(filename = paste0(out_path,"LRB_density_umap.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(t1_1)
dev.off()

null device 
          1 

tiff(filename = paste0(out_path,"ABC_density_umap.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")

print(t2_1)
dev.off()

null device 
          1 

if (!file.exists(paste0(out_path,"virus_integrated_data_080425.rds"))) {
  saveRDS(Sobj, file = paste0(out_path,"virus_integrated_data_080425.rds"))
}

sessionInfo()

R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Los_Angeles
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] scales_1.4.0           ggpointdensity_0.2.0   future_1.58.0          Azimuth_0.5.0         
 [5] shinyBS_0.61.1         patchwork_1.3.1        tidyr_1.3.1            stringr_1.5.1         
 [9] msigdbr_25.1.0         enrichplot_1.26.6      clusterProfiler_4.14.6 harmony_1.2.3         
[13] Rcpp_1.1.0             dplyr_1.1.4            ggplot2_3.5.2          Seurat_5.3.0          
[17] SeuratObject_5.1.0     sp_2.2-0              

loaded via a namespace (and not attached):
  [1] fs_1.6.6                          ProtGenerics_1.38.0              
  [3] matrixStats_1.5.0                 spatstat.sparse_3.1-0            
  [5] bitops_1.0-9                      DirichletMultinomial_1.48.0      
  [7] TFBSTools_1.44.0                  httr_1.4.7                       
  [9] RColorBrewer_1.1-3                tools_4.4.2                      
 [11] sctransform_0.4.2                 R6_2.6.1                         
 [13] DT_0.33                           lazyeval_0.2.2                   
 [15] uwot_0.2.3                        rhdf5filters_1.18.1              
 [17] withr_3.0.2                       gridExtra_2.3                    
 [19] progressr_0.15.1                  cli_3.6.5                        
 [21] Biobase_2.66.0                    spatstat.explore_3.4-3           
 [23] fastDummies_1.7.5                 EnsDb.Hsapiens.v86_2.99.0        
 [25] shinyjs_2.1.0                     labeling_0.4.3                   
 [27] sass_0.4.10                       spatstat.data_3.1-6              
 [29] readr_2.1.5                       ggridges_0.5.6                   
 [31] pbapply_1.7-2                     Rsamtools_2.22.0                 
 [33] yulab.utils_0.2.0                 gson_0.1.0                       
 [35] DOSE_4.0.1                        R.utils_2.13.0                   
 [37] parallelly_1.45.0                 BSgenome_1.74.0                  
 [39] rstudioapi_0.14                   RSQLite_2.4.1                    
 [41] generics_0.1.4                    gridGraphics_0.5-1               
 [43] BiocIO_1.16.0                     gtools_3.9.5                     
 [45] ica_1.0-3                         spatstat.random_3.4-1            
 [47] googlesheets4_1.1.1               GO.db_3.20.0                     
 [49] Matrix_1.7-3                      S4Vectors_0.44.0                 
 [51] abind_1.4-8                       R.methodsS3_1.8.2                
 [53] lifecycle_1.0.4                   yaml_2.3.10                      
 [55] SummarizedExperiment_1.36.0       rhdf5_2.50.2                     
 [57] qvalue_2.38.0                     SparseArray_1.6.2                
 [59] Rtsne_0.17                        grid_4.4.2                       
 [61] blob_1.2.4                        promises_1.3.3                   
 [63] shinydashboard_0.7.3              pwalign_1.2.0                    
 [65] crayon_1.5.3                      miniUI_0.1.2                     
 [67] ggtangle_0.0.7                    lattice_0.20-45                  
 [69] cowplot_1.2.0                     annotate_1.84.0                  
 [71] GenomicFeatures_1.58.0            KEGGREST_1.46.0                  
 [73] pillar_1.11.0                     knitr_1.50                       
 [75] fgsea_1.32.4                      GenomicRanges_1.58.0             
 [77] rjson_0.2.23                      future.apply_1.20.0              
 [79] codetools_0.2-19                  fastmatch_1.1-6                  
 [81] glue_1.8.0                        ggfun_0.1.9                      
 [83] spatstat.univar_3.1-3             data.table_1.17.6                
 [85] vctrs_0.6.5                       png_0.1-8                        
 [87] treeio_1.30.0                     spam_2.11-1                      
 [89] cellranger_1.1.0                  poweRlaw_1.0.0                   
 [91] gtable_0.3.6                      assertthat_0.2.1                 
 [93] cachem_1.1.0                      xfun_0.52                        
 [95] Signac_1.14.0                     S4Arrays_1.6.0                   
 [97] mime_0.13                         survival_3.5-0                   
 [99] gargle_1.5.2                      RcppRoll_0.3.1                   
[101] fitdistrplus_1.2-4                ROCR_1.0-11                      
[103] nlme_3.1-162                      ggtree_3.14.0                    
[105] bit64_4.6.0-1                     RcppAnnoy_0.0.22                 
[107] GenomeInfoDb_1.42.3               bslib_0.9.0                      
[109] irlba_2.3.5.1                     KernSmooth_2.23-20               
[111] seqLogo_1.72.0                    SeuratDisk_0.0.0.9021            
[113] colorspace_2.1-1                  BiocGenerics_0.52.0              
[115] DBI_1.2.3                         tidyselect_1.2.1                 
[117] bit_4.6.0                         compiler_4.4.2                   
[119] curl_6.4.0                        hdf5r_1.3.12                     
[121] DelayedArray_0.32.0               plotly_4.11.0                    
[123] rtracklayer_1.66.0                caTools_1.18.3                   
[125] lmtest_0.9-40                     rappdirs_0.3.3                   
[127] digest_0.6.37                     goftest_1.2-3                    
[129] presto_1.0.0                      spatstat.utils_3.1-4             
[131] rmarkdown_2.29                    RhpcBLASctl_0.23-42              
[133] XVector_0.46.0                    htmltools_0.5.8.1                
[135] pkgconfig_2.0.3                   MatrixGenerics_1.18.1            
[137] fastmap_1.2.0                     ensembldb_2.30.0                 
[139] rlang_1.1.6                       htmlwidgets_1.6.4                
[141] UCSC.utils_1.2.0                  shiny_1.11.1                     
[143] farver_2.1.2                      jquerylib_0.1.4                  
[145] zoo_1.8-14                        jsonlite_2.0.0                   
[147] BiocParallel_1.40.2               GOSemSim_2.32.0                  
[149] R.oo_1.27.1                       RCurl_1.98-1.17                  
[151] magrittr_2.0.3                    GenomeInfoDbData_1.2.13          
[153] ggplotify_0.1.2                   dotCall64_1.2                    
[155] Rhdf5lib_1.28.0                   ape_5.8-1                        
[157] babelgene_22.9                    reticulate_1.42.0                
[159] stringi_1.8.7                     zlibbioc_1.52.0                  
[161] MASS_7.3-58.2                     plyr_1.8.9                       
[163] parallel_4.4.2                    listenv_0.9.1                    
[165] ggrepel_0.9.6                     CNEr_1.42.0                      
[167] deldir_2.0-4                      Biostrings_2.74.1                
[169] splines_4.4.2                     tensor_1.5.1                     
[171] hms_1.1.3                         BSgenome.Hsapiens.UCSC.hg38_1.4.5
[173] igraph_2.1.4                      spatstat.geom_3.4-1              
[175] RcppHNSW_0.6.0                    reshape2_1.4.4                   
[177] stats4_4.4.2                      TFMPvalue_0.0.9                  
[179] XML_3.99-0.18                     evaluate_1.0.4                   
[181] tzdb_0.5.0                        JASPAR2020_0.99.10               
[183] httpuv_1.6.16                     RANN_2.6.2                       
[185] purrr_1.0.4                       polyclip_1.10-7                  
[187] SeuratData_0.2.2.9002             scattermore_1.2                  
[189] xtable_1.8-4                      restfulr_0.0.16                  
[191] AnnotationFilter_1.30.0           RSpectra_0.16-2                  
[193] tidytree_0.4.6                    later_1.4.2                      
[195] googledrive_2.1.1                 viridisLite_0.4.2                
[197] tibble_3.3.0                      aplot_0.2.8                      
[199] memoise_2.0.1                     AnnotationDbi_1.68.0             
[201] GenomicAlignments_1.42.0          IRanges_2.40.1                   
[203] cluster_2.1.4                     globals_0.18.0